/*==================================================================
FILE NAME: Figure_7.do
CREATED: 1 July 2025
====================================================================*/

**Figure 7: Heterogeneity in the Likelihood of a Complaint Investigation and Notice of Violation following a Complaint


/* Set directory if working independently through code
if c(username)=="" { //insert username
	global rootdir "" // insert root path
	global processed_data "$rootdir/processed_data" 
	global figures "$rootdir/output/figures"
	global point_estimates "$rootdir/output/point_estimates"
	// Define global paths for replication package
} 
*/
set scheme modern

use "$processed_data/facility_characteristics_clean.dta", clear
duplicates drop RN ZipCode, force
keep RN ZipCode

merge m:1 ZipCode using "$processed_data/zipcode_demographics_new.dta"
keep if _m == 3
drop _m

gen RN_id = substr(RN,3,.)
destring RN_id, replace
drop RN
save "$processed_data/RN_zipcodes.dta", replace


use "$processed_data/Air_Panel.dta", clear

drop if never_air_inv == 1

* Create time variable
egen t = group(year month) 
xtset RN_id t

* Generate difference variables
forv h = 0/12 {
	gen p_air_complaint_inv_`h' = f`h'.p_air_complaint_inv - l1.p_air_complaint_inv
}

forv h = 2/12 {
	gen p_air_complaint_inv_neg`h' = l`h'.p_air_complaint_inv - l1.p_air_complaint_inv
}

forv h = 0/12 {
	gen p_air_nov_`h' = f`h'.p_air_nov - l1.p_air_nov
}

forv h = 2/12 {
	gen p_air_nov_neg`h' = l`h'.p_air_nov - l1.p_air_nov
}


egen RN_year = group(RN year)

cap drop b u d se Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0

gen b_urban = 0
gen se_urban = 0
gen u_urban = 0
gen d_urban= 0

merge m:1 RN_id using "$processed_data/RN_zipcodes.dta"
keep if _m == 3
drop _m

** URBAN **

gen above_median_urban = 0
sum urban_share, d
replace above_median_urban = 1 if urban_share > `r(p50)'
replace above_median_urban = . if urban_share == .

gen p_air_incident_urban = p_air_incident*above_median_urban 

foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_complaint_inv_`h'  p_air_incident p_air_incident_urban if urban_share !=., absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'

lincom p_air_incident + p_air_incident_urban

replace b_urban = r(estimate) if Years == `h'
replace se_urban = r(se) if Years == `h'
replace u_urban = r(estimate) + 1.96*r(se) if Years == `h'
replace d_urban = r(estimate) - 1.96*r(se) if Years == `h'

}

foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_complaint_inv_neg`h'  p_air_incident p_air_incident_urban if urban_share !=., absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'

lincom p_air_incident + p_air_incident_urban
replace b_urban = r(estimate) if Years == -`h'
replace se_urban = r(se) if Years == -`h'
replace u_urban = r(estimate) + 1.96*r(se) if Years == -`h'
replace d_urban = r(estimate) - 1.96*r(se) if Years == -`h'


}
preserve
keep if Years != .
keep b u d se Years Zero b_urban se_urban u_urban d_urban
export delimited "$point_estimates/Figure_7_Panel_E.csv", replace
graph set window fontface "Times New Roman"
twoway ///
    (rarea u_urban d_urban Years, color(gs2) fintensity(20) lwidth(0) lpattern(solid)) ///
    (rarea u d Years, color(gs10) fintensity(20) lwidth(0) lpattern(solid)) ///
    (line b_urban Years, lcolor(gs3) lpattern(solid) lwidth(medium)) ///
    (line b Years, lcolor(gs6) lpattern(dash) lwidth(medium)) ///
    (line Zero Years, lcolor(gs8) lpattern(solid)) ///
    , ///
    xlabel(-12(1)12, nogrid labsize(vlarge)) ///
    ylabel(0(0.1)0.6, labsize(vlarge)) ///
    legend(order(3 "Above Median Urban Share" 4 "Below Median Urban Share") ///
           pos(6) col(2) size(vlarge)) ///
    ytitle("{&Delta} P(Complaint Investigation)", size(vlarge)) ///
    xtitle("Month", size(vlarge)) ///
    graphregion(color(white)) ///
    plotregion(color(white)) ///
    xsize(8.6)

graph export "$figures/Figure_7_Panel_E.pdf", replace
restore

** Above Median White **

cap drop b u d se Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0

gen b_am_white = 0
gen se_am_white = 0
gen u_am_white = 0
gen d_am_white= 0

gen above_median_white = 0
sum white_share, d
replace above_median_white = 1 if white_share > `r(p50)'
replace above_median_white = . if white_share == .

gen p_air_incident_whiteam = p_air_incident*above_median_white 

foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_complaint_inv_`h'  p_air_incident p_air_incident_whiteam if white_share != ., absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'

lincom p_air_incident + p_air_incident_whiteam
replace b_am_white = r(estimate) if Years == `h'
replace se_am_white = r(se) if Years == `h'
replace u_am_white = r(estimate) + 1.96*r(se) if Years == `h'
replace d_am_white = r(estimate) - 1.96*r(se) if Years == `h'

}

foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_complaint_inv_neg`h'  p_air_incident p_air_incident_whiteam if white_share != ., absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'

lincom p_air_incident + p_air_incident_whiteam
replace b_am_white = r(estimate) if Years == -`h'
replace se_am_white = r(se) if Years == -`h'
replace u_am_white = r(estimate) + 1.96*r(se) if Years == -`h'
replace d_am_white = r(estimate) - 1.96*r(se) if Years == -`h'


}
preserve
keep if Years != .
keep b u d se Years Zero b_am_white se_am_white u_am_white d_am_white
export delimited "$point_estimates/Figure_7_Panel_C.csv", replace
graph set window fontface "Times New Roman"
twoway ///
    (rarea u_am_white d_am_white Years, color(gs2) fintensity(20) lwidth(0) lpattern(solid)) ///
    (rarea u d Years, color(gs10) fintensity(20) lwidth(0) lpattern(solid)) ///
    (line b_am_white Years, lcolor(gs3) lpattern(solid) lwidth(medium)) ///
    (line b Years, lcolor(gs6) lpattern(dash) lwidth(medium)) ///
    (line Zero Years, lcolor(gs8) lpattern(solid)) ///
    , ///
    xlabel(-12(1)12, nogrid labsize(vlarge)) ///
    ylabel(0(0.1)0.6, labsize(vlarge)) ///
    legend(order(3 "Above Median White Share" 4 "Below Median White Share") ///
           pos(6) col(2) size(vlarge)) ///
    ytitle("{&Delta} P(Complaint Investigation)", size(vlarge)) ///
    xtitle("Month", size(vlarge)) ///
    graphregion(color(white)) ///
    plotregion(color(white)) ///
    xsize(8.6)

graph export "$figures/Figure_7_Panel_C.pdf", replace
restore


** Above Median Poverty Rate **

sum adult_pov_rate, d
gen above_med_pov = 0
replace above_med_pov = 1 if adult_pov_rate > `r(p50)'
replace above_med_pov = . if adult_pov_rate == .

gen p_air_incident_pov = p_air_incident*above_med_pov 

cap drop b u d se Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0

gen b_pov = 0
gen se_pov = 0
gen u_pov = 0
gen d_pov = 0


foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_complaint_inv_`h'  p_air_incident p_air_incident_pov if adult_pov_rate != ., absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'

lincom p_air_incident + p_air_incident_pov 
replace b_pov = r(estimate) if Years == `h'
replace se_pov = r(se) if Years == `h'
replace u_pov = r(estimate) + 1.96*r(se) if Years == `h'
replace d_pov = r(estimate) - 1.96*r(se) if Years == `h'

}

foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_complaint_inv_neg`h'  p_air_incident p_air_incident_pov if adult_pov_rate != ., absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'

lincom p_air_incident + p_air_incident_pov 
replace b_pov = r(estimate) if Years == -`h'
replace se_pov = r(se) if Years == -`h'
replace u_pov = r(estimate) + 1.96*r(se) if Years == -`h'
replace d_pov = r(estimate) - 1.96*r(se) if Years == -`h'


}
preserve
keep if Years != .
keep b u d se Years Zero b_pov se_pov u_pov d_pov
export delimited "$point_estimates/Figure_7_Panel_A.csv", replace
graph set window fontface "Times New Roman"
twoway ///
    (rarea u_pov d_pov Years, color(gs2) fintensity(20) lwidth(0) lpattern(solid)) ///
    (rarea u d Years, color(gs10) fintensity(20) lwidth(0) lpattern(solid)) ///
    (line b_pov Years, lcolor(gs3) lpattern(solid) lwidth(medium)) ///
    (line b Years, lcolor(gs6) lpattern(dash) lwidth(medium)) ///
    (line Zero Years, lcolor(gs8) lpattern(solid)) ///
    , ///
    xlabel(-12(1)12, nogrid labsize(vlarge)) ///
    ylabel(0(0.1)0.6, labsize(vlarge)) ///
    legend(order(3 "Above Median Poverty Rate" 4 "Below Median Poverty Rate") pos(6) col(2) size(vlarge)) ///
    ytitle("{&Delta} P(Complaint Investigation)", size(vlarge)) ///
    xtitle("Month", size(vlarge)) ///
    graphregion(color(white)) ///
    plotregion(color(white)) ///
    xsize(8.6)

graph export "$figures/Figure_7_Panel_A.pdf", replace

restore

** NOV **

* Urban *

cap drop b u d se Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0

cap drop b_urban se_urban u_urban d_urban 

gen b_urban = 0
gen se_urban = 0
gen u_urban = 0
gen d_urban= 0

foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_nov_`h'  p_air_incident p_air_incident_urban if urban_share!=., absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'

lincom p_air_incident + p_air_incident_urban
replace b_urban = r(estimate) if Years == `h'
replace se_urban = r(se) if Years == `h'
replace u_urban = r(estimate) + 1.96*r(se) if Years == `h'
replace d_urban = r(estimate) - 1.96*r(se) if Years == `h'

}

foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_nov_neg`h'  p_air_incident p_air_incident_urban if urban_share!=., absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'

lincom p_air_incident + p_air_incident_urban
replace b_urban = r(estimate) if Years == -`h'
replace se_urban = r(se) if Years == -`h'
replace u_urban = r(estimate) + 1.96*r(se) if Years == -`h'
replace d_urban = r(estimate) - 1.96*r(se) if Years == -`h'


}
preserve
keep if Years != .
keep b u d se Years Zero b_urban se_urban u_urban d_urban
export delimited "$point_estimates/Figure_7_Panel_F.csv", replace
graph set window fontface "Times New Roman"
twoway ///
    (rarea u_urban d_urban Years, color(gs2) fintensity(20) lwidth(0) lpattern(solid)) ///
    (rarea u d Years, color(gs10) fintensity(20) lwidth(0) lpattern(solid)) ///
    (line b_urban Years, lcolor(gs3) lpattern(solid) lwidth(medium)) ///
    (line b Years, lcolor(gs6) lpattern(dash) lwidth(medium)) ///
    (line Zero Years, lcolor(gs8) lpattern(solid)) ///
    , ///
    xlabel(-12(1)12, nogrid labsize(vlarge)) ///
    ylabel(-0.02(0.02)0.1, labsize(vlarge)) ///
    legend(order(3 "Above Median Urban Share" 4 "Below Median Urban Share") ///
           pos(6) col(2) size(vlarge)) ///
    ytitle("{&Delta} P(Notice of Violation)", size(vlarge)) ///
    xtitle("Month", size(vlarge)) ///
    graphregion(color(white)) ///
    plotregion(color(white)) ///
    xsize(8.6)

graph export "$figures/Figure_7_Panel_F.pdf", replace
restore


cap drop b u d se Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0

cap drop b_am_white se_am_white u_am_white d_am_white

gen b_am_white = 0
gen se_am_white = 0
gen u_am_white = 0
gen d_am_white= 0


** Above Median White **

foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_nov_`h'  p_air_incident p_air_incident_whiteam white_share if white_share != ., absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'

lincom p_air_incident + p_air_incident_whiteam
replace b_am_white = r(estimate) if Years == `h'
replace se_am_white = r(se) if Years == `h'
replace u_am_white = r(estimate) + 1.96*r(se) if Years == `h'
replace d_am_white = r(estimate) - 1.96*r(se) if Years == `h'

}

foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_nov_neg`h'  p_air_incident p_air_incident_whiteam if white_share != ., absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'

lincom p_air_incident + p_air_incident_whiteam
replace b_am_white = r(estimate) if Years == -`h'
replace se_am_white = r(se) if Years == -`h'
replace u_am_white = r(estimate) + 1.96*r(se) if Years == -`h'
replace d_am_white = r(estimate) - 1.96*r(se) if Years == -`h'


}
preserve
keep if Years != .
keep b u d se Years Zero b_am_white se_am_white u_am_white d_am_white
export delimited "$point_estimates/Figure_7_Panel_D.csv", replace
graph set window fontface "Times New Roman"
twoway ///
    (rarea u_am_white d_am_white Years, color(gs2) fintensity(20) lwidth(0) lpattern(solid)) ///
    (rarea u d Years, color(gs10) fintensity(20) lwidth(0) lpattern(solid)) ///
    (line b_am_white Years, lcolor(gs3) lpattern(solid) lwidth(medium)) ///
    (line b Years, lcolor(gs6) lpattern(dash) lwidth(medium)) ///
    (line Zero Years, lcolor(gs8) lpattern(solid)) ///
    , ///
    xlabel(-12(1)12, nogrid labsize(vlarge)) ///
    ylabel(-0.02(0.02)0.1, labsize(vlarge)) ///
    legend(order(3 "Above Median White Share" 4 "Below Median White Share") pos(6) col(2) size(vlarge)) ///
    ytitle("{&Delta} P(Notice of Violation)", size(vlarge)) ///
    xtitle("Month", size(vlarge)) ///
    graphregion(color(white)) ///
    plotregion(color(white)) ///
    xsize(8.6)
graph export "$figures/Figure_7_Panel_D.pdf", replace
restore


** Above Median Poverty Rate **

cap drop b u d se Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0

cap drop b_pov se_pov u_pov d_pov 

gen b_pov = 0
gen se_pov = 0
gen u_pov = 0
gen d_pov = 0


foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_nov_`h'  p_air_incident p_air_incident_pov if above_med_pov != ., absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'

lincom p_air_incident + p_air_incident_pov 
replace b_pov = r(estimate) if Years == `h'
replace se_pov = r(se) if Years == `h'
replace u_pov = r(estimate) + 1.96*r(se) if Years == `h'
replace d_pov = r(estimate) - 1.96*r(se) if Years == `h'

}

foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_nov_neg`h'  p_air_incident p_air_incident_pov if above_med_pov != ., absorb(RN_year t) cluster(RN_id)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'

lincom p_air_incident + p_air_incident_pov 
replace b_pov = r(estimate) if Years == -`h'
replace se_pov = r(se) if Years == -`h'
replace u_pov = r(estimate) + 1.96*r(se) if Years == -`h'
replace d_pov = r(estimate) - 1.96*r(se) if Years == -`h'


}
preserve
keep if Years != .

export delimited "$point_estimates/Figure_7_Panel_B.csv", replace
graph set window fontface "Times New Roman"
twoway ///
    (rarea u_pov d_pov Years, color(gs2) fintensity(20) lwidth(0) lpattern(solid)) ///
    (rarea u d Years, color(gs10) fintensity(20) lwidth(0) lpattern(solid)) ///
    (line b_pov Years, lcolor(gs3) lpattern(solid) lwidth(medium)) ///
    (line b Years, lcolor(gs6) lpattern(dash) lwidth(medium)) ///
    (line Zero Years, lcolor(gs8) lpattern(solid)) ///
    , ///
    xlabel(-12(1)12, nogrid labsize(vlarge)) ///
    ylabel(-0.02(0.02)0.1, labsize(vlarge)) ///
    legend(order(3 "Above Median Poverty Rate" 4 "Below Median Poverty Rate") ///
           pos(6) col(2) size(vlarge)) ///
    ytitle("{&Delta} P(Notice of Violation)", size(vlarge)) ///
    xtitle("Month", size(vlarge)) ///
    graphregion(color(white)) ///
    plotregion(color(white)) ///
    xsize(8.6)
graph export "$figures/Figure_7_Panel_B.pdf", replace
restore
